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ABSTRACT 

Aims. By reevaluating a 13-month stretch of Ulysses SWICS H pickup ion measurements near 5 AU close to the ecliptic right after 
the previous solar minimum, this paper presents a determination of the neutral interstellar H density at the solar wind termination 
shock and implications for the density and ionization degree of hydrogen in the LIC. 

Methods. The density of neutral interstellar hydrogen at the termination shock was determined from the local pickup ion production 
rate as obtained close to the cut-off in the distribution function at aphelion of Ulysses. As shown in an analytical treatment for the 
upwind axis and through kinetic modeling of the pickup ion production rate at the observer location, with variations in the ionization 
rate, radiation pressure, and the modeling of the particle behavior, this analysis turns out to be very robust against uncertainties in 
these parameters and the modeling. 

Results. Analysis using current heliospheric parameters yields the H density at the termination shock equal to 0.087 + 0.022 cirT 3 , 
including observational and modeling uncertainties. 

Key words. 



1. Introduction 

Neutral interstellar gas of the local interstellar cloud (LIC) pen- 
etrates into the inner heliosphere as a neutral wind due to the rel- 
ative motion between the Sun and the LIC. Apparently, the Sun 
is found near the boundary of a warm, relatively dilute cloud of 
interstellar gas, possibly with a sig nificant gradient in the ion- 
ization fraction o f H and He (e.g. ICheng & Bruhweilerl Il990t 
IWolffetal.|[i999t ISlavin & Frischl |2002l) within a very struc- 
tured surrounding (e.g. reviews bv lCox & Revnoldsll987tlFriichl 
fl995h . In a companion paper within this special section, Frisch 
and Slavin (2008) lay out how the physical parameters and com- 
position of the LIC at the location of the Sun, as derived from in- 
situ observations and from absorption line measurements, con- 
strain the ionization state and radiation environment of the LIC. 
In situ observations of the two main constituents of the LIC, H 
and He, have been obtained with increasing accuracy, starting 
with t he analysis of backscattered solar Lyman-q intensi ty sky 
maps dBertaux & Blamondll97UlThomas & Krassal[l971l) for H 
as wel l as with rocket-borne dParesce et alJI 1974b and satellite- 
borne dWeller & Meier|[l974l) observations of interstellar He us- 
ing the solar He I 58.4 nm line. The optical diagnostics was fol- 
lowed by discovery of pickup io ns for He (Mobius etaDfl98lh 
and for H (Gloeckler et afl 19921) and finally by direct neutral He 
observations dWitte et al.lll993l) . 

Such in-situ diagnostics, even at 1 AU, is made possible by 
the neutral gas flow deep into the inner heliosphere. Through 
the interplay between this wind, the ionization of the neutrals 
upon their approach to the Sun, and the Sun's gravitational field 
(distinctly modified by radiation pressure for H) a character- 
istic flow pattern and density structure is formed, with a cav- 



ity close to the Sun and gravitational focusing on the down- 
wind side (for all species except H). The basic understanding 
of the related heliosphere — LIC i nterac ti on ha s been summa- 
rized i n early reviews by lAxfordl d!972l) ; iFahrl d 19741) : iHolzerl 
dl977l) : lThomasl dl978l) . While He provides us with almost com- 
pletely unbiased information about the physical parameters of 
the LIC since it enters the heliosphere unimpeded, the abun- 
dance of H and O, along with other species, is significantly de- 
pleted, their speed decreased, and their temperature incr eased 
through charge exchange in the heliospheric in t erface (Fahr 
199lURuciriski et al.lll993l:llzmodenov et al.lll999tlMuller et al 
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2000; Izmodenov et al. A consolidation of the physical 

parameters of interstellar He, including the flow velocity vec- 
tor relative to the Sun, as determined from neutral gas, pickup 
ion, and UV backscattering o bservations, was ac hieved through 
the effort of an ISSI Team (Mobi us et alj |2004, and references 
therein), thus leading to a benchmark for the physical parame- 
ters of the LIC. This paper is part of a follow-up effort within 
an ISSI Team to also consolidate the determination of the LIC H 
density. The determination of the H density in the LIC proper not 
only involves a measurement inside the heliosphere, but is also 
dependent on the filtration of H in the heliospheric boundary. 
Therefore, consolidating the observational results concentrates 
on the determination of the H density at the termination shock, 
which still requires taking into account of the dynamics of the 
flow into the inner heliosphere as well as ionization effects. This 
paper deals with the determination of the H density from the 
pickup ion observations made with Ulysses SWICS. The effort 
to obtain the density from mass-loading of the solar wind by H 
pickup ions and its resulting slow down at large distances from 
the Sun is descr ibed in the paper b y Richard son et al.l d2008l this 
volume), while IPrvor et all d2008l this volume) discuss a deter- 
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mination of the H density based on the reduction of the modu- 
lation of the UV backscatter signal with distance from the Sun. 
To illustrat e the state of the m odel-dependent H density value 
in the LIC iMiiller et al.l d2008l this volume) compare different 
global models of the heliosphere and their results for the dis- 
tances of the key boundary structures and the filtration factor, 
which also connects the inner heliosph ere observations to the 
ionization state of the LIC, discussed by Sl avin & Frischl (f2008, 
this volume). 

In their previous work iGloeckler & Geissl d2001al) used 
pickup ion fluxes as observed at 5 AU with Ulysses SWICS, 
the charge exchange rates from SWOOPS, and a Vasyliunas & 
Siscoe distribution function to deduce the local neutral H den- 
sity; they then used a hot interstellar gas model with the ioniza- 
tion rate significantly modified by electron ionization to deduce 
the density at the termination shock. In the present paper we use 
the same data set, but follow a complementary approach. 

After discussing previous derivations of the local neutral gas 
density and its extrapolation to the termination shock at the be- 
ginning of section 2 we present an alternative approach. We 
make use of the fact that the ionization rate in the pickup ion 
production appears both as production rate of PUI and as loss 
rate of the parent neutral gas population and any variations bal- 
ance close to the aphelion of Ulysses. As a consequence, the PUI 
production rate is almost exactly proportional to the H density at 
the termination shock. In the same section we illustrate this be- 
havior in a simplified analytical model that applies to the upwind 
region. 

In section 3 we simulate the local PUI production rate, start- 
ing with the density in the interstellar medium, compare it with 
the observations, and confirm the robustness of this approach by 
varying the parameters. We start with Monte-Carlo simulations 
of the flow through the heliospheric interface for two different 
LIC parameter sets, which result in two different H densities at 
the nose of the termination shock. In a second step, we hand 
these results over to a 3D time dependent test-particle code to 
calculate the H densities and H + PUI production rates at Ulysses 
during the observation interval, while accurately taking into ac- 
count losses and radiation pressure along the trajectories of in- 
terstellar gas. We find the density at the termination shock and 
the LIC parameters that fit the observations best by interpolat- 
ing between the two initial models. We study the response of 
the resulting PUI production rates to variations in the ionization 
rate, radiation pressure, and details in the modeling. In section 4, 
we present the results and show that our method is very robust 
against uncertainties of these parameters, and, in fact, against 
details of the simulations. 
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Fig. 1. The interface correction: ratio of results of the Moscow 
Monte Carlo model for the geometry of Ulysses H PUI observa- 
tions to the results of a sum of classical hot models (two popula- 
tions), evaluated for identical parameters of the gas at the nose of 
the termination shock and identical ionization rate and radiation 
pressure as used in the MC model. The dotted horizontal line 
marks the heliospheric correction value equal to 1, the shaded 
area corresponds to the range of Ulysses heliocentric distances 
during the observations. 

Ulysses locations during the observations period 
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Fig. 2. Change of Ulysses position during the observations. The 
horizontal axis shows the heliocentric distance, the left-hand ver- 
tical scale heliolatitude, and the right-hand vertical scale ecliptic 
latitude. The times are indicated at the plot. Ecliptic longitude 
varied from 153.6° at the beginning of observations interval to 
157.5° at the end, which corresponds to a change from 78.4° to 
81.8° in heliolongitude. 



2. Derivation of the neutral gas density from pickup 
ion fluxes 

The observed pickup ion flux density at any location r in the he- 
liosphere is directly proportional to the local pickup ion source 
function S (r) taken just below the cut-off speed, i.e. at w — 
v/vsw = 1> where v is the pickup ion speed in the rest frame 
of the solar wind, that with respect to the Sun moves with vsw- 
The source function is given by 



S(r)=# on (r)«(r) 



(1) 



where /3 lon (r) is the total local ionization rate and n (r) the local 
interstellar neutral gas density. Hence determining the local neu- 
tral gas density requires the knowledge of the ionization rate, 
and the uncertainty of the derived density is directly related to 
the uncertainty with which the ionization is known. 



The situation appears even worse in the attempt to derive 
the interstellar gas density at the termination shock from obser- 
vations in the inner heliosphere. For most interstellar species, 
except for He, ionization has already significantly depleted the 
local density at least where good quantitative observations of in- 
terstellar pickup ions have been available so far. In this way the 
local neutral gas density is also dependent on the ionization rate, 
with a depletion that typically scales exponentially with the in- 
verse of the distance from the Sun and could be written for any 
short local stretch as: 



n (r) ~ exp 



Aon (r) r 



(2) 



Only the proportionality is important here for the arguments 
made below and not any constants, such as a, that can be ad- 
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Measured absolute PUI production rates 



8x10"" 
7x 10" 1 ' 




2.5 3 3.5 4 4.5 5 

heliocentric distance 



Fig. 3. Absolute production rates of H PUI as a function of he- 
liocentric distances, obtained from the observed PUI spectrum 
during the interval discussed in the paper. Shown are the rates 
after normalizing the spectrum to the production rate at 1 AU 
and then multiplying by 1/r 2 . 



justed for normalization. In addition, the average ionization rate 
relevant for the depletion of the neutral gas (referred to as the 
loss rate) may be different from that responsible for pickup ion 
generation (production rate) because of the different time scales 
involved, and, in particular, because H is also subject to radi- 
ation pressure and its variations. This combination makes any 
determination of the density of neutral gas at the termination 
shock dependent both on the modeling and on the knowle dge of 
heliospheric parameters (e.g. iRucihski & Bzowskil ll996). This 
is certai nly true for the H PUI observations taken w ith Ulysses 
SWICS (lGloecklerill99rtlGloeckler & Geissll2001bh at or inside 
5.3 AU. 

2.1. Previous neutral gas density determination 

To minimize this influence, Gloeck ler & Geissl d2001bl) used 
an approach that a) relies on a long-term averaging of data at 
Ulysse s, b) uses the PUI transport model by lVasyliunas & Siscod 
dl976l) . and c) simultaneously determines the total ionization 
rate from the slope of the pickup ion velocity distribution, 
which reflects the radial distribution of the neutral gas inside 
the observer distanc e. A similar approach had been taken by 
Mobius et al. (1988) for the determination of the He density 
from observations at 1 AU. By making use of the fact that 
He + pickup ions are solely created by charge exchange with 
solar wind He 2+ and of the Ulysses SWICS capability to si- 
multaneously observe He 2+ pickup ion and solar wind fluxes at 
~ 5 AU, where interst ellar He is not significantly depleted yet, 
IGloeckler et all (fl997l) were able to obtain a He density whose 
uncertainty only depends on the knowledge of the charge ex- 
change cross section and is independent of the absolute calibra- 
tion of the observing instrument. For H at least two ionization 
processes contribute substantially, solar wind charge exchange 
and UV ionization, so that the rate cannot be easily eliminated 
from the analysis by simultaneous measurements. This remain- 
ing uncertainty in the total ionization rate, which translates into 
the uncertainty of the den sity, is exemplified by the fact that 
iGloeckler & G eiss (2001a) had to invoke a rather high electron 
ionization rate of 2.4 x 10~ 7 s _1 at 1 AU (with a distance de- 
pendence that is stronger than 1/r 2 ) to explain the pickup ion 
velocity distribution. 



Recently, IGloeckler et al.l d2008l) determined the densities 
of interstellar H, N, O, Ne, and Ar at the termination shock 
using their abundances relative to He in the energetic tails of 
the ion distributions in the heliosheath, as obtained with both 
Voyager LECP sensors with a relative uncertainty of +10%. To 
arrive at the absolute densities, they used the interstellar He 
density derived from Ulysses SWICS He ++ pickup ion mea- 
surements (see previous paragraph) of 0.015 + 0.002 cm -3 
(IGloeckler & Geissll2004T; IGloeckler et al.ll2004bl) . which is also 
the consensus value for the interstellar He density based on these 
picku p ion, direct neutral gas, and UV backscattering observa- 
tions^ (Mobius et al. 2004). Combining these two observations, 
IGloeckler et al.l d2008l) arrived at an H density at the termina- 
tion shock of . 08 cm ' 3 . With the two uncertainties cited by 
IGloeckler et al.l (120081) combined as independent contributions, 
the resulting uncertainty for the H density is obtained equal to 
+0.013 cirT 3 , or 17%. While this method is insensitive to uncer- 
tainties in the absolute geom etric factor of the Ulys ses SWICS 
instrument of +25%, cited by IGloeckler et al.l (12008), which ap- 
plies to the direct H pickup ion observation method, the deter- 
mination of the abundances may be subject to additional uncer- 
tainties in the production rates for the different species that were 
used to infer the neutral species ratios. 



2.2. Alternative neutral gas density determination, 
minimizing the influence of uncertainties in the 
ionization rate 

In the following discussion we take a different approach, which 
will minimize the influence of uncertainties in the ionization rate 
on the resulting neutral gas density at the termination shock. We 
will illustrate this behavior in a simplified analytical treatment 
in this section, before we show with parameter variations in rig- 
orous simulations, discussed in the next section, that this also 
holds for the actual Ulysses observations and even extends to 
variations in the radiation pressure. 

As becomes obvious from the relation between the neutral 
gas density and the pickup ion source function, a linear de- 
pendence of the resulting neutral gas density on the ionization 
rate (according to Eq.(Q}) remains valid even at the termina- 
tion shock. Conversely, an exponential behavior due to deple- 
tion by ionization (see Eq.([2jl) prevails in the inner heliosphere, 
massively overcompensating the linear dependence of the source 
function on the ionization rate. Between these two extreme lo- 
cations must lie a place where the effects of the ionization rate 
on the source function and thus on the observed pickup ion flux 
cancel. Here, the observed quantity is strictly proportional to the 
neutral density at the termination shock and - to the first order 
- does not depend on the choice of the ionization rate any more, 
if the adopted ionization rate is not too different from the correct 
value. We will explore this behavior below for inflow on the he- 
liospheric upwind axis for which an analytical solution can be 
found. 

In our derivation we make the reasonable assumption that ra- 
diation pressure fully compensates solar gravity: y. — 1. This as- 
sumption is justified during the observation period, as discussed 
in detail in the Appendix. Even a significant reduction of yiz has 
very little influence on the final result, as is shown with the sim- 
ulations in section 3. For // close to 1 the interstellar inflow speed 
remains almost constant as a function of distance from the Sun 
and equal to the interstellar bulk flow speed at the termination 
shock Vtsm- Then the density of the neutral interstellar gas is 
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reduced along the upwind axis according to 
dn (r) 



At 



= -Aon n (r) . 



With a constant inflow speed we can use: 
- dr = Vism dr 



(3) 



(4) 



(with the gradient directed inward). The ionization rate, both so- 
lar wind charge exchange and UV combined, varies as 



AonO) = Aon,E rl/r 2 ' 



(5) 



where Aon,£ is the ionization rate at r# = 1 AU. Combining 
Eq.© through Eq.© together leads to 



dn 0) Aon,£ r\ dr 



n(r) 



Vism 



After logarithmic integration this yields: 
In 

which is equivalent to 

n (r) = hq exp 



w(r) _ Aon£ 4 
«0 Vism?" ' 



VlSM?" J 



A 

«o exp I 



(6) 



(7) 



(8) 



i.e., the density falls off exponentially with the typical penetra- 
tion distance A = Aon,£ r^/VisM, as the gas approaches the Sun. 
Consequently, the pickup ion source function is 



S 0) = Aon,£ Cy\ «o exp 



/^ion,£ Y £■ 



Vism'" 



which now depends linearly on the neutral density no at the ter- 
mination shock and in two ways on the ionization rate Aon,£- 
Equation (O is not dependent on Aon,£ anymore if dS /dAon,£ = 
for a given distance r. This condition yields 



1 = Aon,E r|/ViSM r 



(10) 



or the effects of the ionization rate cancel exactly at the pen- 
etration distance A, i.e. at the edge of the hydrogen cavity in 
the heliosphere which, by definition, is the geometric location 
of the surface where the local density is equal to 1/e of the 
de nsity at "infinity". For the ionization rate of 5.5 10~ 7 s~' used 
by iGloeckler & Geissl d2001iJ) and an interstellar inflow speed 
of 22 km/s and temperature of ~ 12 000 K at the termination 
shock, as results for the combined primary and secondary distri- 
bution s of interstellar H from global modeling (Izmodenov et alJ 
2003 b]) and from SOHO SWAN observations dOuemerais et alJ 
1999tlLallement et al.ll2005llCosta et al.ll 19991) . the point of per- 
fect compensation is at 3.8 AU in the upwind direction, where 
the cavity edge is closest to the Sun. At crosswind, where the 
Ulysses observations were made, the cavity ends at ~ 5.5 AU, 
and in the downwind direction at a still larger distance. 

At these distances from the Sun the local density of neutral 
gas depends linearly on the local ionization rate and on the den- 
sity at the termination shock n®. Consequently, the source func- 
tion of PUI, as defined in Eq.(Q~|i, depends linearly on no- 
Through kinetic simulations that include all important ef- 
fects in the inner heliosphere, such as ionization and ra- 
diation pressure, we will demonstrate in the following sec- 
tion that the Ulysses observat ions in 1997 through 199 8 be- 
tween 5 and 5.4 AU, used by IGloeckler & Geissl d2001al) and 



llzmodenov et alJ d2004l) . were indeed made in a region where 
the effects of a potentially not so well known ionization rate can- 
cel and also that uncertainties in the radiation pressure as well 
as effects from different treatments of the particle behavior in 
the modeling are minimal. Henceforth, we will make use of the 
same data set to derive a refined value for the interstellar H den- 
sity at the TS, which is robust against remaining uncertainties in 
the heliospheric parameters that control the local density distri- 
bution such as ionization rates and radiation pressure. 

3. Simulations and comparison with the 
observations 

In the following we will model the interstellar H distribution at 
the locations of the observations, starting from the pristine in- 
terstellar medium. We will reproduce the observed production 
rate of H pickup ions while taking into account all known and 
relevant heliospheric processes and their current uncertainties as 
much as possible. 

Calculation of the local H density at the point of Ulysses 
observations involves length scales that span two orders of mag- 
nitude, from the scale of the typical penetration distance of H 
A =i 3 AU, to the size of the heliosphere > 100 AU. Likewise, 
the time scales have a comparable range - from weeks for solar 
UV illumination structures and rotation period to ~ 40 years of 
the H travel time from the pristine LIC to the observation point. 
To cover such ranges in a single simulation at a sufficient level 
of detail is beyond the reach of current computer resources avail- 
able. 



(9) 3.1. Description of numerical models 



Therefore, separate simulations were performed on two different 
spatial scales and levels of detail, and then combined. In a first 
step, large scale configuration of the heliosphere was established 
for two sets of parameters for the LIC and solar output aver- 
aged over a large time interval, using t he Moscow Monte Carlo 
(MC) code dBaranov & Malamal fl993). In a next step, modifi- 
cations of the neutral interstellar hydrogen flow due to its inter- 
action with the solar environment inside the termination shock 
were simulated using the Warsaw test particle 3D and time de- 
pendent kinetic code dRuciriski & Bzowski 1995; Bzowsk i et al.l 
119971 l2002l:lTarnopolski & Bzowskill2007l) . 

The model adopted in the MC simulation is static and axi- 
ally symmetric, assuming a constant and spherically symmetric 
solar particle and radiation output. The MC model is used to in- 
fer the physical parameters of H at the nose of the termination 
shock, i.e. the flow that will finally reach the inner solar sys- 
tem, where the observations are made. As extensive ly discussed 
in the past dBaranov et al.|[l998t Ilzmodenovll2000l) . neutral in- 
terstellar gas at the termination shock can be approximated to 
some extent by two populations (primary and secondary), fea- 
turing Maxwellian distributions shifted in velocity by specific 
bulk flow values. The primary population represents the original 
interstellar neutral H gas, while the secondary population is cre- 
ated between the heliopause and heliospheric bow shock due to 
charge exchange between interstellar atoms and the compressed 
and heated plasma that flows around the heliopause. These two 
populations are taken as boundary conditions at the termination 
shock for the test-particle model of the inner heliosphere. 

The current Warsaw test-particle 3D and time-dependent 
model requires a known distribution function far away from the 
Sun, which is invariable in time and homogeneous in space. 
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Consequently, the two neutral-gas populations are Maxwellians 
characterized by temperature, bulk velocity, and density. 

As demonstrated by iBaranov et al.l d!998l) : llzmodenovl 
(2001); Izmod enov et al.l d2001l) . such assumptions are a consid- 
erable simplification. In reality, neither the distribution functions 
of the two populations are Maxwellian, nor their macroscopic 
parameters are hom ogeneous in space. On the other hand, it was 
also demonstrated (Izm odenov et alj|2005l) that 30% sinusoidal 
variations in the spherically symmetric solar wind density dur- 
ing the solar activity cycle induce only ~ 10% variations in the 
thermal populations of neutral H at the termination shock. Since 
simulations including more realistic variations of the solar wind 
evolution are unavailable, and given tremendous computer bur- 
den of the MC time-dependent model, these variations were as- 
sumed to be negligible for this study. 

Deviations of the parameters of the distribution function 
from homogeneity at the termination shock result in systematic 
differences in the local densities returned by the Moscow MC 
and Warsaw test-particle models. To assess the robustness of our 
approach to such deviations, we calculated first the H density at 
the geometric location of the PUI data with the use of the MC 
model and then, adopting identical values of the radiation pres- 
sure and ionization rates and the parameters of the primary and 
secondary populations at the termination shock, with the Warsaw 
test-particle model in its axially symmetric, static mode. The ra- 
tio of the density profiles returned by the two models is shown 
in FigQ] While the two results significantly deviate from each 
other closer to the Sun along the Ulysses PUI accumulation line, 
the MC and test particle models return almost identical results 
between 3.5 and 5 AU. Generally, the hot model overestimates 
the density by a factor which increases towards the Sun. Since in 
the present study we are interested only in the agreement of the 
models at Ulysses location (~ 5 AU from the Sun), where the 
results of the MC and test particle models agreed to 1 %, we feel 
justified to adopt in the analysis the values for the PUI produc- 
tion rates at Ulysses calculated with the use of the Warsaw test 
particle model. 

Earlier versions of the model were d escribe d by 

Rucihski & Bzowsl dl d!995l) : iBzowski et all d 19971) and 
Bzo wski et al.l d2002l). The present v e rsion of the model was 
described by Tarnopolski & Bzowski (2007) and includes the 
following effects: 

- the ionization rate, composed of photoionization, charge- 
exchange and electron-impact ionization rates; the net ion- 
ization rate varies with heliolatitude and its radial profile dif- 
fers only slightly from 1/r 2 ; 

- the radiation pressure can include the dependence of its mag- 
nitude on the radial velocity of individual atoms with respect 
to the Sun; the net intensity is also variable with time, which 
results in non-Keplerian trajectories of the atoms; 

- the inflowing neutral gas is split into two populations of 
thermal atoms with different parameters at the termination 
shock. 

We constrained the simulations by aligning the relevant 
model parameters with available data wherever possible. Only 
in lack of available data proxies and models were used to infer 
the necessary parameters. In order to not disrupt the flow of the 
discussion, details of the radiation pressure and ionization rate 
models are presented in the Appendix. In the following section, 
we will introduce the use of the pickup ion data and how the 
computation has been adapted to them. 



3.2. Pickup ion data and appropriate computation mesh 

For our comparison we use the local product ion rate of H + 
picku p ions as measured by SWICS/Ulysses (Gloeckler et al. 
119921) from 1997.285 until 1998.310, when Ulysses was cross- 
ing the ecliptic plane at aphelion of its orbit, going from ~ 4.95 
to ~ 5.4 AU and descending from 20° to 0° ecliptic latitude, 
whic h corresponds to an interv al from +12° to -6° heliolatitude 
("c.f. iGloeckler & Ge iss 2001b). The spectrum that was used to 
derive the production rate as a function of distance from the 
Sun, a quantity that is directly related to the pickup ion flux 
and thus very close to the observable in this measurement, (as 
presented in Fig. [3} is an ensemble average of many individ- 
ual spectra registered during the time interval mentioned. The 
individual spectra were selected so that the phase space den- 
sity in the suprathermal tails (v/v sw > 2.4) was minimum. This 
eliminated contributions from shocks in the solar wind. The so- 
lar wind proton peak was corrected for the instrument dead- 
time effects. The averaged spectrum was fitted by forward mod- 
eling using the cl assical hot model and the theory of pickup 
ion transport from Vasvliunas & Siscoe] d 19761) . The best fit was 
obtained for the following set of the hot model parameters: 
= 0.9, p = 6.1 x 10~ 7 s -1 , v TS = 22 km/s, T JS = 12000 K, 
«ts =0.1 cirT 3 . This procedure returned a spectrum that then 
was expressed as the PUI production rate normalized to 1 AU 
from the Sun. The absolute PUI production rates as a function of 
distance from the Sun that are shown in Fig. [3] were obtained by 
multiplying this fitted spectrum by by 1/r 2 . The PUI production 
rate used to derive the H density at the TS was taken from the 
portion 1.92 < v/v sw < 2.01 in the spectrum, which corresponds 
to the distances marked in Fig.Q] Its magnitude was determined 
to be equal to 7.26 x 10~ 1() cirT 3 s _! with an experimental un- 
certainty of +25%, almost entirely attributed to the systematic 
uncer tainty in the geometric factor of SWICS (IGloeckler et al.l 
2008), because statistical fluctuations are negligible for the long 
time average used here. 

Fluctuations in the pickup ion fluxes due to transport effects 
on shorter time scales do not affect the pickup ion production rate 
based on the full observation interval. In the following we use 
this value to determine the density of neutral interstellar H at the 
termination shock in a comparison with the simulated pickup ion 
production rate, while varying the ionization rate and radiation 
pressure within the range of recent observations. 

In order to adapt the simulations to the long observation in- 
terval with changing locations of Ulysses we performed the cal- 
culations of the local hydrogen density on a mesh of = 16 
points distributed evenly in time along the Ulysses orbit during 
the observation interval. The densities of the two populations 
were computed separately for each of the two populations at the 
16 points and then combined to obtain the net local density for a 
given moment of time: 

»; (r) = n P rU ( r > &h *i) + "sec,; (r, A/, </>;, f,) , (11) 

where « pr j,,-, « S ec,/ are local densities of the primary and secondary 
populations, r is the Ulysses radial distance, f,- is the i-th time 
moment and A;, 4>i are Ulysses heliocentric coordinates at f, (see 
Fig.|2]i. Further, the production rates calculated at each of the 16 
points were averaged. Therefore, the resulting mean production 
rate for a given simulation was calculated as follows: 

N 

ySprod = Yj ["'■ (r ' Ai > fc« fi) ^ (r ' fi) J /N > (12) 

!=1 

where /3 is the net local ionization rate of neutral hydrogen. 
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3.3. Calculations 

For the boundary conditions in the Local Interstellar Cloud 
(LIC) we adopted the gas inflow direction, bulk velocity, and gas 
temperature as deriv ed by the ISSI Wor king Group on Neutral 
Interstellar Helium (Mobius et al. 2004) from in situ observa- 
tions of neutral interstellar H e atoms (Witt e 2004), from mea- 
surements of He pickup ions dGloeckler et al. l2004al). and from 
UV o bservations of the heliospheric He glow ( Vallerga et al. 
|2004|) . The upwind direction adopted in the simulations was 
A B = 254.68°, <t> B = 5.31° in the B 1950.0 ecliptic coordinates; 
bulk velocity vg = 26.4 km s _1 ; and temperature Tb = 6400 K. 
The density of neutral He in the LIC was adopted as equal to 
0.015 cm" 3 , and based on the He ionization degree in the LIC, 
inferred by I Wolff et ail d 19991) on the level of ~ 30 - 40%, the 
density of He + in the LIC was taken equal to 0.008 cm 3 . 

In our simulations we used two paramete r sets for LIC H 
(referre d to by 1 . and 2. below), ad opted from Izmo denov et al.l 
d2003bl) and llzmodenov et al.l (l2003al) . which include: 

1 . LIC proton density n p = 0.06 crrT 3 and neutral gas density 
«e =0.18 cm 3 , which yields the H ionization degree in the 
LIC equal to 25%; the contribution of He + to the net plasma 
density in the LIC would hence be on the level of 12%. 

2. LIC proton density n p = 0.032 cirT 3 and neutral gas density 
«h = 0.2 cirT 3 , with the same density of He + , which yields 
the H ionization degree of only 14%. 

Running the Moscow MC model as the first step of the 
simulation resulted in the following parameters of the primary 
and secondary populations at the nose of the termination shock, 
which were adopted in the second step of the simulations with 
the use of the Warsaw test-particle code: 

1. Primary: WTS.prf = 0.1925 n H = 0.03465 cm 4 , Vxs.pri 
1 .08 v B = 28.512 km sec" 1 , r TS ,pri = 6020 K. 
Secondary: nxssec = 0.3345 «h = 0.06021 cirT 3 , VTS.sec = 
0.71vb = 18.744 km sec" 1 , r TS ,sec = 16300 K. 

The resulting net density at the termination shock was thus 
equal to «h,ts,i = 0.53 n H = 0.095 cirT 3 . 

2. Primary: «TS P ri = 0.2926 n H = 0.05852 cm 4 , Vxs.pri 
1 .07 v B = 28.248 km/s; r T s,p ri = 6100 K; 

Secondary: nxssec = 0.2934 n H = 0.05868 cirT 3 , v T s se c 
0.7 v B = 18.48 km/s; r TS , S ec = 16500 K. 
The resulting net density at the termination shock was thus 
equal to nH.TS.n = 0.59 «h = 0.117 cm 4 and the contribu- 
tion of He + to the plasma density in the LIC would be much 
higher, namely 20%. 

In the second step of the simulations, the H densities and H + 
PUI production rates at Ulysses were calculated with the Warsaw 
test-particle model using Eqs ( fTTT i and dTZI ), with the parameters 
of the two populations at the TS as boundary conditions. In gen- 
eral, these models can be subdivided into two groups, depending 
on the treatment of temporal variations in radiation pressure and 
ionization rate. In one of the groups, the parameters were cal- 
culated as monthly averages for the 16 time intervals f,, and the 
test-particle program was run in the static mode for each time 
interval f;, with the parameters pertinent to f; ("instantaneous" 
simulation). In the second group, the parameters were taken as 
"smooth" analytic models, as presented in the Appendix, and 
the test-particle program was run in its time-dependent mode for 
each time interval f,- ("smooth" simulation). 

With such an approach we could assess the influence of 
short-time fluctuations on the averaged result, which turned out 



to be negligible (though differences between "smooth" and "in- 
stantaneous" values at f,'s were on the order of 30%). The re- 
sult obtained in a comparison of monthly averages and a smooth 
temporal variation on timescales of one year and longer also jus- 
tifies our simplification to ignore completely time scales shorter 
than one month. Since the "instantaneous" simulation was much 
less computer-intensive, most of the further tests were run in the 
"instantaneous" mode. 

Additional variable elements in this part of the simulations 
were the treatment of the radiation pressure (either dependent on 
or independent of the radial velocity of the atoms), the inclusion 
or exclusion of a latitudinal anisotropy in the ionization rate, and 
inclusion or exclusion of the electron impact ionization. 
Variation of heliospheric control parameters 

In a next step of the simulations, the robustness of the H + 
PUI production rate at Ulysses against variations in the absolute 
values and treatment of the radiation pressure, ionization rate, 
and in the modeling strategies was tested. We have repeated the 
second step of the simulations with various ionization rate and 
radiation pressure values within the limits of their observational 
uncertainties. 

These tests included reducing the ionization rate compared 
with the most recent compilation. Since the charge exchange 
rates as recorded on Ulysses close to the ecliptic (normalized 
to 1 AU) appear systematically lower than the values obtained 
from the OMNI 2 time series, we repeated step 2 simulations 
with the equatorial charge exchange values yS eqtl (see Appendix) 
reduced to Ulysses measurements, which resulted in an overall 
reduction of the ionization rate by ~ 25%. 

Since the absolute calibration of the solar Lyman-o- input has 
changed appreciably since the beginning of the observations in 
the 1970-ties, when it was believed that the effective radiation 
pressure factor a t solar minimum was n 0.7 and at solar maxi- 
mum// 1 (e.g. IVidal-Madiarll 1 9751: iTobiska et aill 1 997l) . com- 
pared to present -day values ^ 1 at solar minimum and ^ 1 .5 at 
solar maximum (Woods et al. 2000; TobiskaetaJJ&OOO), we re- 
peated step-2 calculations for the line-integrated flux reduced by 
factors 0.9, 0.8, and 0.7, as illustrated in Fig]4] We also checked 
for robustness against uncertainties in the solar Lyman-o- line 
profile by performing simulations with a simplified flat line- 
center (hence, the radiation pressure independent of atom radial 
radial velocities) and a self-reversed profile (with radiation pres- 
sure dependent on the radial velocity). 

As can be seen from the previous discussion, the variations 
in the ionization rate and radiation pressure introduced into our 
simulations were not chosen arbitrarily, and they did not exceed 
uncertainties related to changes in calibrations and increased 
measurement sophistication. 

To also assess the sensitivity of the results to uncer- 
tainties in the cross section for charge exchange, we re- 
peated the simulations f or the first LIC parameter set with the 
Maher & Tinsl ey| (| 19771) cro ss section replac ed with that from 
Lindsay & Stebbingsld20()5t) . As discussed bv lFahr et all d2007l) . 
the two formulae agree to a few percent in the supersonic so- 
lar wind regime, but differ up to ~ 40% for low collision en- 
ergies, pertinent to the region between the heliopause and the 
bow shock, where the primary interstellar population loses a por- 
tion of its atoms and the secondary population is created due to 
charge exchange with protons. As illustrated by iBaranov et al.l 
d 19981) . the change in coupling between the protons and H atoms 
results in different proportions between the primary and sec- 
ondary populations at the termination shock. Our simulations 
showed that apart from the changes in the individual densities of 
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Fig. 4. Model/data ratios S,,i/S obs and S,, 2 /Sobs of the H + PUI 
production rates at Ulysses as a function of the reduction factor 
of the solar radiation pressure relative to the currently adopted 
value. The group of lines with filled symbols (1) corresponds 
to the simulations performed with the density at the termination 
shock equal to 0.095 cm 3 , the open symbols (2) correspond to 
0.117 cm 4 . Diamonds correspond to the nominal values of the 
ionization rate and radiation pressure dependent on v r , stars to 
the reduced values of the ionization rate and radiation pressure 
dependent on v r , and squares and triangles to radiation pressure 
independent on v r , with the ionization rate, respectively, nomi- 
nal and reduced. The observed PUI production rate was equal to 
7.26 x 10- 10 cnr 3 s" 1 . 



the populations at the termination shock, their remaining param- 
eters, i.e. bulk velocities and temperatures, change very little: 
Primary: «xspri = 0.02592 cm 4 , vjspri = 28.776 km sec -1 , 
r T s, P ri = 5900K, 

Secondary: «TSsec = 0.07020 cm -3 , VTSsec = 18.520 km sec -1 , 
7rs, S ec= 16500 K; 

the net density at the termination shock was equal to 0.096 cm -3 , 
i.e. practically identical as the in the simulation with the old 
cross section formula. 

Hence, the propagation and losses of the two populations 
during their travel from the termination shock to the inner he- 
liosphere were almost identical in the cases of Maher & Tinsley 
and Lindsay & Stebbings formulae, but the input values to the 
simulation inside the termination shock were different. 



4. Results 

The results of the extended simulations are compiled in Fig0]as 
ratios of the simulated production rates to the measured value. 
The simulations show that indeed, as postulated in Section 2, the 
model H + PUI production rate at the location of Ulysses during 
the observation interval is only weakly dependent on the radi- 
ation pressure, ionization rate, and details of modeling of the 
gas density in the inner heliosphere. A change in the ionization 
rate by ~ 25% (between diamonds and stars or triangles and 
squares in Fig|4| results in a change in the density at the ter- 
mination shock of only ~ 2.5%. The combined variation in the 
H + PUI production rate due to details of the ionization rate and 
the modeling approach for radiation pressure does not amount 
to more than 4%. Varying the level of solar Lyman-a' output by 
30% translates into a somewhat larger variation in the production 
rate, on the level of 10%, but again this is substantially weaker 
than the variation in the input. Overall, the amplitude of the vari- 
ations in each of the two input factors is reduced by a factor of 



10 for the ionization rate and 3 for the radiation pressure in the 
resulting variations in the inferred interstellar H density at the 
termination shock. 

This result is also robust against other details of the simula- 
tions, including the solar Lyman-a' line shape (inverted profile 
vs flat), functional form of the charge exchange cross section 
dMaher & Tinslevld 19771) vs lLindsav & Stebbings! d2005l) ). treat- 
ment of temporal variations in radiation pressure and the ioniza- 
tion rate (fully time-dependent vs static with solar-rotation av- 
eraging), presence or absence of the latitudinal variations in the 
ionization rate, and presence or absence of the electron ioniza- 
tion. All these factors affect the production rate at Ulysses during 
the observation period only by a few percent and are not able to 
significantly change its value. 

To calculate the density at the termination shock, we take the 
PUI production rates for the nominal jj values and the /i val- 
ues reduced by 10%, obtained from the first simulation (filled 
symbols in Fig.|4]i, and we we calculate the density at the termi- 
nation shock «h,ts from the formula: 



«H,TS - 



«H,TS,1 
(S i,\/S bs) 



(13) 



where the angular brackets denote arithmetic mean, S 0DS is the 
PUI production rate from the Ulysses observations discussed in 
this paper, S i.i are the production rates from the first (1) set 
of simulations, and «h,ts,i is the density at TS assumed in the 
first simulation (denoted with index 1). The «h,ts value obtained 
from this calculation is equal to 0.089 + 0.003 cm 4 . 

The robustness of this result is supported by the derivation 
of the density at the TS from the other simulation (2), for which 
the TS density for the nominal set of ionization and radiation 
pressure values was assumed to be 0.117 cm 4 (the results are 
shown as open symbols in Fig. @). Deriving the TS density 
analogous to the description in the previous paragraph yields 
»h,ts = 0.086 + 0.003 cm 4 , i.e. very close to the previous value 
although the starting TS density in the simulation was higher by 
35%. 

Hence, by taking the average of the two values, we arrive at 
the density of interstellar hydrogen at the TS of 0.087 cm 4 + 
25%, where the uncertainty is almost completely dominated by 
the instrumental uncertainty of the absolute geometric factor. 

5. Discussion and conclusions 

We have used the accumulation of the H + pickup ion production 
rate from SWICS/Ulysses over a ~ 13 month period in 1997 - 
1998 at the Ulysses passage through the solar equator plane to 
infer the interstellar H density at the termination shock. By ex- 
tensive simulations we demonstrated that the H + PUI production 
rate in this location of the heliosphere is only weakly dependent 
on the values of solar radiation pressure and neutral H ioniza- 
tion rate, but sensitively depend on the density at the termination 
shock. We have found that the H density at the termination shock 
inferred from these pickup ion production rates is very robust 
against any variations in the ionization rate, radiation pressure, 
and the actual modeling approaches for the density distribution 
in the inner heliosphere. 

In the present analysis we have, for the first time, included 
explicitly both the observational uncertainty and the modeling 
uncertainties. While in our approach the modeling uncertain- 
ties are minimized to a few percent, a larger uncertainty is 
incurred for the observation because absolute flux values are 
used, which results in an uncertainty of the obtained termina- 
tion shock density equal to +25%. In the previous approach 



8 



Bzowski et al.: Density of neutral H from Ulysses pickup ion observations 



bv lGloeckler & Geissl (1200 lab the observational uncertainty was 
minimized by making use of the ratio of the pick up ion and solar 
wind flu x, but the ~ 10% uncertainty quoted in Izmo denov et alj 
(2003b) does not include a range of values for the ionization rate 
and radiation pressure and the uncertainty of the geometric fac- 
tor of the instrument. But, as demonstrated in the previous sec- 
tions, the density at the termination shock scales linearly with the 
observed pickup ion production rate, which is directly related to 
the pickup ion flux and/or distribution function. Hence, any ob- 
servational unc ertainties will transfer line arly into the resulting 
densities. Since ITzmodenov et alj (l2003bl) started from the local 
neutral gas density at Ulysses, any uncertainty in the ionization 
rate will appear approximately linearly in the extrapolated den- 
sity at the termination shock. 

The determination of the PUI production rate at Ulysses near 
the aphelion, on which our derivation of «h,ts is based, is not 
entirely model-free. Although in the present determination of the 
PUI production rate at Ulysses a simple hot model was used for 
forward-modeling of the PUI distribution function, our method 
is robust against simplifications inherent to that kind of modeling 
because it uses a quantity (i.e., the production rate of PUI at 
Ulysses) which weakly depends on details of such modeling. 

The determinat i on of the H density at the termination shock 
bv lGloeckler et all d2008l) . equal to 0.080 + 0.008 cm 3 , is free 
of the uncertainty of the geometric factor of the instrument, but 
is subject to a combination of uncertainties in the determination 
of the He density and that of the He abundance relative to H 
from the Voyager LECP observations, including uncertainties in 
the ratios of the production rates of these species. Nevertheless, 
after including of all uncertainti es, all three app r oaches (i.e. 
the present one and th ose from [Gloeckler & Geiss (2001aj) and 
iGloeckler et alj d2008l) ) should be read with a similar uncertainty 
band. The density value presented here agrees very well with 
the new determination bv lGloeckler et al.l d2008l) . Although our 
value still agrees with the previous determin ation of the H den- 
sity fr om SWICS pickup ion observations bv lGloeckler & Geissl 
(2001a) within their mutual uncertainty bands, the combination 
of the two new results suggest a somewhat lower density than 
0.1cm- 3 . 



Our results also agree comfortably with the TS density val- 
ues found from th e analysis of the heliospheric Lyman-a glow 
dPrvor et al.ll2008l this volum e) and from the solar wind slow- 
down ( Richardson et al.l 120081 this volume) within the uncer- 
tainty bands. It should be noted here that the coupling be- 
tween the neutral and ionized component of the interstellar 
medium between the bow shock and the heliopause appears to 
be somewhat stronger than suggested previously. This is a con- 
sequence of an updated relation for the energy dependence of 
the charge exchange cross-s ection between protons and H atoms 
dLindsav & Steb bings 200J). 



Lyman-cr flux 



The parameters of the interstellar gas in front of the helio- 
spheric bow shock, as a ssumed in simulatio n (1), also seem to 
be robust, as shown by Miilleret al. (2008, this volume), who 
discussed the present status of the modeling of heliospheric in- 
terface and showed that differences in the filtration rate returned 
by different models of the heliosphere evaluated with identical 
initial parameters are about 15% and this result can be adopted 
as the uncertainty of the H density in the CHISM. 
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Fig.A.l. Daily values of line-integrated solar Lyman-cc flux 
(small dots) from SOLAR 2000, Carrington-averaged values 
(thick dots) and fitted models (smooth lines). The upper (bro- 
ken) line corresponds to the approximation of the Carrington- 
averaged time series by Eq. dA.U . The lower (solid) line indicates 
the effective radiation pressure after rescaling to the central band 
with Eq. dA.2| i. The horizontal line corresponds to// = 1. 



Appendix A: Inner-heliospheric environment during 
observations 

A. 1 . Radiation pressure 

Radiation pressure varies with heliocentric distance, radial ve- 
locity of individual atoms, and time. Both of the cases which will 
be discussed below are based on the line- and disk-integrated so- 
lar Lyman-a flux I tot at 1 AU (expressed in photons crrT 2 s -1 ), 
based on daily valu es obtained from the SOLAR 2000 model 
(Tobiskaetai][2000), which are shown for the observation inter- 



val as small dots in Fig. [AT] 

Because it takes years for the interstellar flow to pass even 
only through the inner heliosphere, we ignored variations on 
time scales shorter than one solar rotation period (i.e., a month), 
which is further justified by the fact that two different treatments 
on longer time scales lead to almost identical results. To assess 
the influence of slower temporal variations on the modeling we 
constructed two models: "instantaneous" and "smooth". In the 
"instantaneous" model the appropriate monthly values, shown 
as thick dots in Fig. IA.1I were taken within a time-independent 
model for the entire month. In the " smooth" model, th e variation 
in 7 tot was approximated, following BzowskJ d2001bl) . by the re- 
lation: 

hot (0 = hot,o + ^ (fli cos o>it + bi sin w,f) (A. 1) 
i=i 

with N M = 7. The relevant frequencies and their ampli- 
tudes wer e obtained from the analysis of the Lomb peri- 
odogram (Press & Rybicki 1989) of a time series composed 
of the Carrington-averages of 7 t ot for the time period 1948- 
2004 (results obtained using daily values are almost identical). 
Periodicities shorter than one year and amplitudes smaller than 
0.025 of the strongest harmonic were ignored. The parameters 
obtained are listed in Table lA.fl and the resulting curve is shown 
as the upper broken line in Fig. I A. 1 1 

The solar Lyman-o- line has a two-peaked, self-reversed pro- 
file (Fig. IA.2t . The radiation pressure acting on individual H 
atoms depends on the Doppler shift resulting from their radial 
motion relative to the Sun. Therefore, the Lyman-or spectral flux 
was converted into the radiation pressure factor fi, using either a 
"flat" or a "Doppler" model. 
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Solar Lyman-o- line profiles 



Radiation pressure vs total flux 
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Fig. A.2. Model sol ar Lyman-*? line profiles, based on data from 
Lemai re et al. (2002), for the start and end of the PUI observa- 
tion interval, expressed in the yU units. Thick lines indicate the 
approximate radial velocity range for the H atoms at Ulysses, 
adopted as the bulk velocity +3x the thermal speed. The range 
is used for the "flat" model in Eq. flA.2| >. 



Since most of the H atoms at ~ 5 AU do not exceed 
~ 30 km/s, one can approximate the spectral flux responsi- 
ble for the radiation pressure by averaging the line profile over 
~ +30 km/s about the line center for a few solar line profile 
data sets and relating these values to the routinely -measured 7 tot . 
Based on observations by Lemaire et al. (2002), the averaged 
central Lyman-a flux correlates linearly with the total flux 7 tot , 
as illustrated in Fig. IA.3I and a relation for the "flat" radiation 
pressure factor /ig at can be given as: 



3.473 x 10" I2 / tO f 



0.287. 



(A.2) 



A similar fit was published bv lEmerich et al.l d2005l) for the cen- 
ter of the line profile, i.e. radial velocity of km/s. The "flat" 
approximation for the radiation pressure has been used in con- 
nection both with the "instantaneous" and "smooth" time series 
of the Lyman-a fluxes. 

The "Doppler" model of the radiation pressure is based on 
a f unctional fit to the 9 solar Lyman-a line profiles observed 
by Lemaire et al.l (120021). The model, discussed in detail by 
Tarnopolski & Bzo wskil d2007l) . uses the functional form: 

H (v V , / to t) = A ( 1 + B /tot) exp (-Cv?) (A.3) 
x [l + D exp (Fv r - Gvfj + H exp {-Pv r - Qvfjj 

where v r is the radial velocity of a H atom in km/s. The param- 
eters of the fit are compiled in Table IA.2I Sample fits for the 
beginning and the end of the pickup ion observation interval are 
shown in Fig. IA.2I 

In summary, four baseline models were compared, i.e. the 
"flat" approximation and the "Doppler" model of the radiation 
pressure in combination with the "instantaneous" and "smooth" 
temporal dependences. The absolute variation in the radiation 
pressure discussed in the main text was introduced by multiply- 
ing each model by scaling factors ranging from 0.7 to 1 . 

With these in hand, one can construct the radiation pressure 
for a given time t and radial velocity v, by inserting of the / to t (?) 
value to Eq.( IA.3b , obtained from either the smooth or instanta- 
neous model. 

Hence, a total of 4 baseline models of radiation pressure 
were exercised: two "flat" models, with / to t either instantaneous 
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Fig. A.3. Relation between the line-integrated flux of the solar 
Lyman-ff output, expressed in the //-units, and the coefficient /i 
of radiation pressure acting on interstellar hydrogen atoms in 
the inner heliosphere. Dots correspond to the values obtained 
from averaging of t he 9 solar Lyman-a profiles observed by 
iLemaire et alJ (120021) over ±30 km/s about the line center. The 
line is a fit defined in Eq.( IA.2l ). 



Table A.l. Parameters of Eq. dA.U needed to calculate the line- 
integrated solar Lyman-ff flux for time t in years; / tot ,o = 4.6034x 
10 11 ctrT 2 s- 1 . 



i 


COi 


Clj 


bi 


1 


0.14406 


-1.7673 x 10 1U 


-1.5657 x 10 1U 


2 


0.31182 


-6.3068 x 10" 


1.0732 x 10 10 


3 


0.43745 


-2.1161 x 10" 


8.3570 x 10 9 


4 


0.58427 


5.0496 x 10 10 


7.9792 x 10 10 


5 


0.74027 


-2.0954 x 10 10 


-9.6306 x 10 9 


6 


1.14430 


-1.8691 x 10 10 


-1.3267 x 10 9 


7 


1.96351 


-6.2468 x 10 9 


-3.8329 x 10 9 



Table A.2. Parameters of the model of radiation pressure depen- 
dence on radial velocity v r and total flux / to t defined in Eq.( lA.31 l. 



A = 

D 
H 



2.4543 x 10~ 9 , 
: 0.73879, 
= 0.47817, 



B 
F 
P 



4.5694 x 10~ 4 , 
: 4.0396 x 10~ 2 , 
4.6841 x 10~ 2 , 



C : 
G : 

Q 



3.8312 x 10\ 
: 3.5135 x 10- 4 , 
: 3.3373 x 10" 4 



or smooth, and two "Doppler" models, also with / to t either in- 
stantaneous or smooth. The reduction of radiation pressure dis- 
cussed in the main text was executed by multiplying relevant 
model by scaling factors ranging from 0.9 to 0.7. 



A.2. Ionization processes 

Ionization of H in interplanetary space occurs through a com- 
bination of three processes. Charge exchange with solar wind 
and photoionization by solar EUV are usually two major con- 
tributors and both scale almost perfectly with the inverse square 
of the distance r from the Sun. Electron impact ionization is 
mostly a small contribution that increases closer to the Sun (in- 
side ~ 3 AU), but becomes completely negligible at large dis- 
tances. 
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A.2.1. Charge exchange 

The charge exchange rate within the supersonic solar wind was 
calculated in the standard way (e.g. iBzows kl l2001ah . with the 
use of the widely adopted formula: 



Charge exchange rate 



jScx = 0" cx Osw) «pVsw 



(A.4) 



where vsw is the solar wind speed, «sw solar wind density, and 
cr cx is the reaction cross section, which depends on the relative 
velocity of the particles. At a given heliolatitude the solar wind 
speed was assumed to be constant up to the termination shock, 
with the density to drop off as 1/r 2 . Hence the charge exchange 
rate decreases with 1/r 2 . 

Similarly to our treatment of radiation pressure, two ap- 
proaches were used to model the time variations: "smooth" and 
"instantaneous". As input we used the equatorial rate at 1 AU, 
denoted yS eqtl -. Another variable in the model was the presence 
or absence of latitudinal anisotropy. This resulted in 4 baseline 
models: spherically symmetric "smooth" or "instantaneous" and 
latitude-dependent "smooth" or "instantaneous". In the latitude- 
dependent model the anisotropic part of the model was calcu- 
lated as shown below, for any given time t . In the "instantaneous" 
case the latitudinal structure was "frozen" as for this time, while 
in the "smooth" case a fully time-dependent model was real- 
ized, with the latitudinal rate evolving along the trajectory of the 
atoms. 

Evolution of the equatorial rate 

The approach to model the equatorial rate of charge ex- 
change was similar to the modeling of the line-integrated flux 
of the solar Lyman-a radiation in Section Al, now based on the 
daily values o f the solar wind speed and density from the OMNI- 
2 collection dKing & Papitashvihl 120051) . normalized to 1 AU. 
The charge exchange rate was calculated according to Eq.( lA.4l > 
and subsequently averaged over Carrington periods. The results 
were taken as the "instantaneous" model for the relevant 16 time 
intervals. For the "smooth" model a periodogram analysis was 
performed, which returned a formula similar to Eq. flA.ll ), with 
N C x - 5 periodicities and the remaining parameters collected in 
Table lA.3l The results are shown in Fig. IA.4l 

A comparison of the daily charge exchange rates based on 
OMNI-2 with observations on Ulysses (scaled to 1 AU) showed 
systematic differences, with the Ulysses values lower by ~ 25%. 
This resulted in another set of models, with the rates reduced by 
25% relative to the OMNI-2 values. 
General formula for the 3D charge exchange rate 

To perform a realistic simulation of the Ulysses PUI obser- 
vations one has to include the evolution in the charge exchange 
rate with heliolatitude. Based on observations. lBzows ki (2001a) 
proposed an analytical phenomenological formula for the charge 
exchange rate at 1 AU, repeated here for reader's convenience in 
a slightly modified form: 



yScxOM = (j8pol+fy0) + (/?eqtr(O-A>ol) 



(A.5) 



x exp 



ln J 2(»-fe (t)-(f>s(t) 
\ to (»Ws(0 



(p is the heliographic latitude, n a shape factor which was adopted 
as n — 2, y8p i the mean charge exchange rate at the poles. The 

term (j3 po i + 6p 4>) describes the north-south asymmetry of the po- 
lar rates, and the term exp [ — In 2 ( ) ] provides the lati- 
tudinal evolution. 

Evolution of latitudinal anisotropy 
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Fig. A.4. Comparison of the model equatorial charge exchange 
rates with the Ulysses observations obtained during the PUI 
measurementsinterval. The broken line is the "instantaneous" 
model, the solid line is the "smooth" model (Eq. dA.lt ), and dots 
are the OMNI-2 daily rates. 



Ulysses in situ observations dMcComas et al.1 1 19991) were 
not sufficient to infer the evolution in the latitudinal anisotropy 
and had to be supplemented by observations of the heliospheric 
Lyman-a glow, which is sensitive to the st ructure of the sola r 
wind. As demonstrated by SWAN/SOHO dBertaux et alJI 19991) . 
the glow features a darkening in an ecliptic band during low 
solar activity, nicknamed the heliospheric groove. The groove 
is due to the latitudinal anisotropy of the sola r wind, makin g 
it a good trac e r for the latitudinal structure dBzows kH 120031) . 
Bzo wski et al.l d2003l) exploited this to infer the equatorial-to- 
pole contrast of the charge exchange rates and the latitudinal 
boundaries of the polar regions of reduced rates for selected 
dates between solar minimum and maximum. 

With this information and continuous coverage of the equa- 
torial charge exchange rate, "snapshot pictures" of the charge 
exchange ionization field in the inner heliosphere were worked 
out. 

Gaps in the coverage of the Lyma n-a images were f i lled b y 
observations of polar holes reported bv lHarvev & ReceTvl (|2002), 
who provide a time series of the latitudinal boundaries of polar 
holes between two consecutive solar maxima. 

We took advantage of a ne w linear correlation betw een the 
areas of the polar holes from lHarvev & Recelvl d2002l) S C h,N, 
5 c h,s, defined as follows: 

I cos(0) d(f> dA = 2n (1 - sin (<f> c hjf)) 

T ch,N JO 

cos(0)d0d/l = 27r(l +sin(0 ch , 5 )), (A.6) 

n/2 JO 

an d the areas of r e duced charge exchange rate Sn, Ss> inferred 
by iBzowski et al.l (|2003) from observations of the heliospheric 
glow. (f> c h,N, 4>ch,s are the latitudes of the north and south hole 
boundaries. Both boundaries are shown in Fig |A.5l Since, as in- 
ferred from the observations, these areas vanish at solar max- 
imum, coefficients were fitted separately to the northern and 
southern areas: 



S n - a^S c h,rv 
S s - as S c h,s 



(A.7) 



with a^/ = 3.91 and as - 4.20. 

Sn, Ss are computed from Eqs. dA.7b using the observed po- 
lar hole areas based on the boundaries cp C h,N, <p c h,s ■ Finally, the 
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Fig. A.5. Latitudinal bound aries of coronal holes as reported 
by lHarvev & ReceTvl d2002l) (top and bottom lines) and bound- 
aries of the "equator ial bulge" of the charge exchange rate from 
iBzowski etail (120031) (solid and dotted lines). 

Table A.3. Parameters to calculate the equatorial charge ex- 
change rate for time t in years; /8 e qti-,o = 5.2978 x 10~ 7 s _1 . 



i 


a>j 


a, 


b, 


1 


0.15671 


-6.4582 x 10~ 8 


-4.0790 x 10- 8 


2 


0.39108 


3.3963 x 10- 8 


3.3016 x 10- 9 


3 


0.59997 


-3.8191 x 10- 8 


-1.7056 x 10~ 8 


4 


0.96353 


-2.5358 x 10- 8 


3.9630 x 10-" 


5 


1.21876 


-3.6585 x 10- 8 


-1.9796 x 10- 8 



Table A.4. Parameters of Eq. dA.9t needed to calculate the 
boundaries of polar holes for time t expressed in decimal years 



N/S 


CJ0 


00 


0i 


N 


0.58251 


31.2 


-24.0 


S 


0.58226 


-38.7 


21.5 



boundaries of the reduced ionization rate regions in the northern 
and southern hemispheres emerge from Eq. (IA.8b as follows: 

4>n = arcsin(l - 2Sn) 

cf>s = -arcsin(l - 2S S ) (A.8) 

As shown in Fig |A.6l this result compares well with the 
boundaries inferred from SWAN observations. To simplify com- 
putations, the boundaries of the coronal holes were further ap- 
proximated by: 

4>n,s (t) = 0o + 0i exp [- cos 3 (w 1)\ (A.9) 

with the parameters listed in Table |A. 4"! 

North - south asymmetry of the charge exchange rate 

The north-south asymmetry in the polar charge exchange 
rates was discover ed during the first Fast La t itude Scan by 
Ulyss es in 1995 dMcComas et alj 1 19991 l2000l IBzowski et at] 
2003). Its long-term evolution is still unknown. As an ad hoc so- 
lution, we include this effect with the term {3 po \+6p (p m Eq. dA.5l l, 
where /3 poi = 2.5405 x 10~ 7 and S p = -1.7863 x 10- 10 . 

A.2.2. Photoionization 

Over time scales longer than one solar rotation and at distances 
outside a few dozen solar radii, the photoionization field can 
be treated as spherically symmetric, with the intensity falling 
off with the square of the distance. Departures from spherical 



1990 1992 1994 1996 1998 2000 2002 
time lyr] 

Fig.A.6. Boundaries of the "equatorial bulge" in the charge 
exchange rate based on coronal hole boundaries and Eq.( IA.8l > 
(soli d and dotted l ines), compared with the boundaries inferred 
by IBzowski et al.l (|2003) from the heliospheric Lyman-o' glow 
(broken lines). 
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Fig. A.7. EUV ionization rates: "instantaneous" (thick dots) and 
"smooth" model (continuous line). Small dots represent daily 
values as computed from the 10.7 cm proxy. 

symmetry are about ~ 10% dAuchere et al.l 120051) and were ne- 
glected. 

The photoionization rate /3 p h was treated similarly to the net 
Lyman-o' flux 7 tot , with the "instantaneous" and "smooth" mod- 
els. Measurements of the photoionization rates of H over the so- 
lar cycle are not readily available, so proxies have to be used. A 
reasonable proxy for the H photoionization rate is the solar ra- 
dio flux in the 10.7 cm band (Bzowski 2001b). Thus daily values 
of the Ottawa solar 10 .7 cm flux were u sed, converted to daily 
photoionization rates dBzowskil l2001bl) . and then Carrington- 
averaged, as shown for the "instantaneous" and "smooth" mod- 
eling in Fig. IA.7I The rate relevant for the Ulysses PUI obser- 
vations was about 0.8 x 10~ 7 s _1 , slowly increasing with time. 

When this research was already pretty much advanced, a 
new version of SOLAR 2000 was introduced, with the H pho- 
toionization rates deduced from the proxies-inferred solar spec- 
tra. A comparison of the rates derived in the two ways showed 
that a difference in the amplitudes exists, while the solar cycle- 
averages agree. The ratio of amplitudes is ~ 1 .4. Using the two 
predictions for the simulations of the Ulysses observations, we 
find an agreement of the results within ~ 10%. Because the sim- 
ulations are time consuming and the influence of photoioniza- 
tion is weak, we decided to keep the model based on the 10.7 cm 
proxy. Since the proxies used in the SOLAR 2000 model are 
much more elaborate than the simple 10.7 cm proxy, the H pho- 
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toionization rate based on SOLAR 2000 is probably an improve- 
ment and thus should be used in future studies, in particular 
during solar minimum and maximum, when the differences are 
largest. 



A.2.3. Ionization by electron impact 

Throughout the simulation we used only one model of the 
electron-impact ionization rate. The electron impact rate was 
calculated based on the radial ev olution in the solar wind 
electron temperature as derived by Mars chet al.l (Q989) from 
Helios da t a and on the mean solar wind density taken from 
Koh nleinl d!996l) . The electron density was adopted assuming 
quasi-neutrality of the solar wind, as a sum of the proton density 
and twice the alpha density. 

The electron distribution function can be approximated as a 
bi-Maxwellian, with a warm core and hot halo populations, plus 
an oc casional strahl po pulation along the local magnetic field 
line dPilipp et al.lll987aTlbl) . The contribution of the halo popula- 
tion to the net rate is on the level of a few pe rcent and is an in- 
creasi ng function of the heliocentric distance (Maksimovic et al. 
120051) . The calculations tha t we performed based on the ACE 
electron data presented by McM ullin et al.l d2004t) show that at 
1 AU the ionization rate due to the core population of the so- 
lar wind electrons is equal to about 0.4 x 10~ 7 s and to the 
halo population to less than 0.04 x 10~ 7 s _1 . Our assessment 
showed also that the amplitude of fluctuations in the electron 
ionization rate may reach an order of magnitude, which is much 
more than the long-time variations related to variations over the 
solar cycle. On th e other hand, the electron data from Wind 
dSalem et alj|2003l) lead to an in-ecliptic solar minimum (1995) 
rate of ~ 0.68 x 10~ 7 s _1 and a solar maximum (2000) rate of 
~ 0.73 x 10~ 7 s _1 . Thus assuming a constant rate over the solar 
cycle is a good approximation. 

Because of the lack of electron distribution data during the 
observation interval, the rate was calculated assuming a mono- 
Maxwellian elec t ron d istribution, following the approach by 
Rucihski & Fahri dl99ll) . Since the analytical formula for the 
electron ionization rate is still quite complex and consider- 
ably slows down the simulation, we used an approximate phe- 
nomenological formula, where the heliocentric distance r is ex- 
pressed in AU and the rate (r) at r is expressed in s^ 1 by: 



2 , , / 10.95(lnr- 124.1)(lnr + 6.108) \ 

r Pei(r) = exp 

H w F \ (In r- 7.491)(65.25 + In r(lnr+ 15.63))/ 

(A.10) 

The formula, shown in Fig |A.81 is valid to the distance of ~ 
10 AU, beyond which the electrons are cooled so much that their 
ionization capability is negligible. 

Observations do ne with Ulysses (Phillips et al.l 119951; 
llssautier et al.l Il998h suggest that the electron ionization rate 
is a 3D, time depend ent function of the solar cycle phase. 
McMullin etal. (2004) came to a similar conclusion for the elec- 
tron impact rate of helium. However, the PUI measurements with 
Ulysses were collected outside ~ 2 AU, where the relative con- 
tribution of electron ionization is already very small. Therefore, 
the simplified model presented above is well justified. A model 
of the electron-impact ionization rate relevant for lower helio- 
centric distances, including its latitudinal evolution durin g the 
solar cycle, has recently been presented bv lBzowskH (|2008). 



Electron ionization rate of H 




0.1 0.2 0.5 1 2 5 10 

heliocentric distance (AU) 

Fig. A.8. Radial profile of the electron-impact ionization rate of 
neutral interstellar H used in the simulations (solid line). The 
radial dependence of this rate differs significantly from the 1/r 2 
law valid for charge exchange and photoionization, as illustrated 
by the broken line. 
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Fig.A.9. Net ionization rates scaled to 1 AU as a function of 
time. Dots represent the "instantaneous" model, the solid line is 
the "smooth" model evaluated for 0°heliolatitude. 



A.2.4. Net ionization rate 

The net ionization rates were calculated as a sum of the charge 
exchange, photoionization, and electron impact rates. They are 
presented in Fig |A.9l as a function of time, normalized to the lat- 
itude of ecliptic plane at 1 AU, and in Fig. |A.10| as a function of 
heliolatitude, both for the "instantaneous" and "smooth" mod- 
els. In the simulations we used (i) a spherically symmetric, in- 
stantaneous model, (ii) a 3D (anisotropic) instantaneous model, 
(3) a spherically symmetric smooth model, and (4) a 3D smooth 
model. 
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Fig. A.10. Net ionization rates scaled to 1 AU as a function of he- 
liolatitude. Dots represent "instantaneous" values and solid lines 
the 3D "smooth" evaluated for the beginning and end of the ob- 
servation interval. 
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